Incorporating New Evidence and Uncertainty

Backwards Conversion, Solving for PSA Parameters and Copula-Based PSA Sampling

source("./manifest.r")

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.1     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'kableExtra'


The following object is masked from 'package:dplyr':

    group_rows



Attaching package: 'Matrix'


The following objects are masked from 'package:tidyr':

    expand, pack, unpack


here() starts at /Users/johngraves/Dropbox/teaching/SMDM-Europe-2023

Loading required package: forecast

Registered S3 method overwritten by 'quantmod':
  method            from
  as.zoo.data.frame zoo 

Registered S3 methods overwritten by 'demography':
  method      from 
  print.lca   e1071
  summary.lca e1071

This is demography 2.0.0.9000 



Attaching package: 'MASS'


The following object is masked from 'package:dplyr':

    select



Attaching package: 'expm'


The following object is masked from 'package:Matrix':

    expm

Learning Objectives

Back-convert an existing transition probability matrix to incorporate new health states, strategies, and evidence.
Solve for probabilisitic sensitivity analysis (PSA) distribution parameters.
Sample correlated PSA distributions using copulas.

1. Backwards Conversion

Lectures 2 and 3 emphasized the importance of the transition rate matrix as the central “hub” of a Markov model.

Goal

What Rarely Happens

Adapting a Model Requires a Transition Rate Matrix

1. Backwards Conversion

If starting a model from scratch, can simply define or estimate rates, and use them to construct the rate matrix.

1. Backwards Conversion

What if we have a model that is already defined in terms of a transition probability matrix?
How can we convert back to a rate matrix to add new health states, transition states, accumulators, etc. as needed?
Also useful if we want to keep everything the same, but change the model time step (e.g., see Chhatwal, Jayasuriya, and Elbasha (2016))

1. Backwards Conversion

We will next show you several ways to work backwards.
Boils down to solving for the continuous “generator matrix” for the observed transition probability matrix.

A Word of Warning

As you’ll see, it’s not always going to work perfectly.
If the original transition probability matrix was defined incorrectly (e.g., no jumpover probabilities), the generator matrix may not exist.
Identifiability is a more general issue, however.

Working Example: HIV Model

Transition probability matrix:

mP = matrix(c(0.721, 0.202, 0.067, 0.010,
                0.000, 0.581, 0.407, 0.012,
                0.000, 0.000, 0.750, 0.250,
                0.000, 0.000, 0.000, 1.000), 
              nrow=4, byrow=TRUE,
              dimnames=list(c("Healthy", "LowCD4", "AIDS", "Death"),
                            c("Healthy", "LowCD4", "AIDS", "Death")))
mP

        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.581 0.407 0.012
AIDS      0.000  0.000 0.750 0.250
Death     0.000  0.000 0.000 1.000

From ‘Decision Modelling for Health Economic Evaluation’ by Briggs, Claxton, Sculpher.

Three Methods to Solve for the Generator

A. Compute the matrix logarithm of the transition probability matrix $P$.
B. Maximum likelihood-based approach.
C. Convert the supplied transition probabilities to rates one-by-one.
(You can always try 2+ of these methods to see how they compare.)

A. Matrix Logarithm of P

In mathematical terms, the generator matrix is the matrix logarithm of the transition probability matrix.
A matrix has a logarithm if and only if if it is invertible.

A. Matrix Logarithm of P

In mathematical terms, the generator matrix is the matrix logarithm of the transition probability matrix.
A matrix has a logarithm if and only if if it is invertible.

\[ R = \log P \]

A. Matrix Logarithm of P

The $\log$ can be found using spectral or eigenvalue decomposition.
If $V$ is a matrix where each column is an eigenvector of $P$, then

A. Matrix Logarithm of P

The $\log$ can be found using spectral or eigenvalue decomposition.
If $V$ is a matrix where each column is an eigenvector of $P$, then

\[A' = V^{-1} A V\]

A. Matrix Logarithm of P

The $\log$ can be found using spectral or eigenvalue decomposition.
If $V$ is a matrix where each column is an eigenvector of $P$, then

\[A' = V^{-1} A V\]

\[\log P = V (\log A') V^{-1}\]

A. Matrix Logarithm of P

V  <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R  <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3)

A. Matrix Logarithm of P

V  <- eigen(mP)$vectors
iV <- solve(V)
Ap <- iV %*% mP %*% V
lAp <- diag(log(diag(Ap)), nrow(Ap), ncol(Ap))
R  <- V %*% lAp %*% iV
R[abs(R) < 1e-6 ] <- 0
dimnames(R) = dimnames(mP)
round(R,3)

        Healthy LowCD4   AIDS  Death
Healthy  -0.327  0.311  0.002  0.013
LowCD4    0.000 -0.543  0.615 -0.072
AIDS      0.000  0.000 -0.288  0.288
Death     0.000  0.000  0.000  0.000

A. Matrix Logarithm of P

Rexp <-
  R %>%
  round(., 3) %>%
  data.frame() %>% 
  rownames_to_column() %>% 
  mutate_all(as.character) %>% 
  mutate(Death = cell_spec(Death,"html",color = ifelse(rowname=="LowCD4","red","black"))) 

Rexp %>% 
  kable(col.names = c("","Healthy","LowCD4","AIDS","Death"),escape = FALSE) %>% 
  kable_styling("hover", "striped", full_width = F)  %>% 
  kable_paper()

	Healthy	LowCD4	AIDS	Death
Healthy	-0.327	0.311	0.002	0.013
LowCD4	0	-0.543	0.615	-0.072
AIDS	0	0	-0.288	0.288
Death	0	0	0	0

Note there is a negative rate from LowCD4 Death!
This is a transition from Death LowCD4 and should be moved to the other side of the matrix.

A. Matrix Logarithm of P

Rexp %>% 
  kable(col.names = c("","Healthy","LowCD4","AIDS","Death"),escape = FALSE) %>% 
  kable_styling("hover", "striped", full_width = F)  %>% 
  kable_paper()

	Healthy	LowCD4	AIDS	Death
Healthy	-0.327	0.311	0.002	0.013
LowCD4	0	-0.543	0.615	-0.072
AIDS	0	0	-0.288	0.288
Death	0	0	0	0

The negative rate implies an identifiability issue.
Note rates from HealthyAIDS, HealthyDeath are relatively small, implying these are from statistical noise in observation.

A. Matrix Logarithm of P

expm::logm Higham (2008) method returns same as eigenvalue method.
Tweaking rates to get a proper model is ad-hoc and difficult.

A. Matrix Logarithm of P

Higham (2008) method returns same as eigenvalue method. - expm::logm()
Tweaking rates to get a proper model is ad-hoc and difficult.

R = expm::logm(mP)

R <- logm(mP)
dimnames(R) = dimnames(mP)
round(R,3)

        Healthy LowCD4   AIDS  Death
Healthy  -0.327  0.311  0.002  0.013
LowCD4    0.000 -0.543  0.615 -0.072
AIDS      0.000  0.000 -0.288  0.288
Death     0.000  0.000  0.000  0.000

B. Multi-State Model Based Approach

Imposition of structural assumptions and fitting to pseudo-data derived from original Markov model can result in reasonable rate estimates.
We could assume that a patient with HIV at this point in time only gets sicker and external causes of death are negligible. Healthy LowCD4 AIDS Death constrains model.
We use the reported data from the original model at year 1.

B. Multi-State Model Based Approach

Steps:

Run the original model for 1+ cycles to obtain the Markov trace.
Construct psuedo-data from the resulting trace based on a cohort with reasonable size (e.g., 1,000 patients)
Estimate transition hazards in the pseudo-data based on a multistate model.¹
Use the estimated transition hazards to construct the rate matrix.

1. Run the original model

tr0 <- # Initial state occupancy.
  c("Healthy" = 1000,"LowCD4" = 0, "AIDS" = 0, "Death" = 0)
tr1 <- # State occupancy after one cycle.
  tr0 %*% mP 

tr <- # Bind together into a data frame. 
    rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    select(cycle,everything())
tr

  cycle Healthy LowCD4 AIDS Death
1     0    1000      0    0     0
2     1     721    202   67    10

2. Construct psuedo-data

Idea: create data for 1,000 simulated “patients” followed for two cycles.
Counts in the Markov trace govern state occupancy in each cycle.
- In first cycle, all 1,000 are in Healthy state.
- In second cycle, 721 remain in Healthy while 202 are in LowCD4, etc.

2. Construct psuedo-data

tr <- 
    rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    gather(state,count,-cycle) %>% 
    mutate(count = round(count,0)) %>% 
    arrange(cycle) %>% 
    lapply(.,rep,.$count) %>% 
    cbind.data.frame() %>% 
    as_tibble() %>% 
    mutate(state = factor(state,
      levels = c("Healthy","LowCD4","AIDS","Death"))) %>% 
    arrange(cycle,state) %>% 
    mutate(ptnum = rep(1:1000,2)) %>% 
    select(-count) %>% 
    arrange(ptnum,cycle) %>% 
    mutate(state = as.numeric(state)) %>% 
    select(ptnum,cycle,state)

tr

# A tibble: 2,000 × 3
   ptnum cycle state
   <int> <dbl> <dbl>
 1     1     0     1
 2     1     1     1
 3     2     0     1
 4     2     1     1
 5     3     0     1
 6     3     1     1
 7     4     0     1
 8     4     1     1
 9     5     0     1
10     5     1     1
# ℹ 1,990 more rows

3. Fit Multistate Model

Next, fit a multistate model to the pseudo-data.
The multistate model will use likelihood-based methods to estimate transition intensities among health states in the data.
These transition intensities are the rates that you can use in the rate matrix!

3. Fit Multistate Model

Define the state table for the multistate model

library(msm)


Attaching package: 'msm'

The following object is masked from 'package:flexsurv':

    qgeneric

statetable.msm(state, ptnum, data=tr)

    to
from   1   2   3   4
   1 721 202  67  10

Check that cell transition counts match the trace.

# Markov trace
rbind.data.frame(tr0, tr1) %>%
    mutate(cycle = c(0, 1)) %>% 
    select(cycle,everything())

  cycle Healthy LowCD4 AIDS Death
1     0    1000      0    0     0
2     1     721    202   67    10

3. Fit Multistate Model

Initial rate guesses are based on the eigenvalue method, i.e., result of logm(mP).

Q.init <- rbind(rbind(c(0, 0.3, 0,   .0001),    
                      c(0, 0,   0.6, 0.01),    
                      c(0, 0,   0,   0.3),  
                      c(0, 0,   0,   0)))
dimnames(Q.init) = dimnames(mP)

Fit the model

hiv.msm <- 
  msm(state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)

3. Fit Multistate Model

hiv.msm


Call:
msm(formula = state ~ cycle, subject = ptnum, data = tr, qmatrix = Q.init)

Maximum likelihood estimates

Transition intensities
                  Baseline                          
Healthy - Healthy -0.32712 ( -3.680e-01, -2.907e-01)
Healthy - LowCD4   0.32702 (  1.354e-01,  7.898e-01)
Healthy - Death    0.00010 (  0.000e+00,        Inf)
LowCD4 - LowCD4   -0.64477 ( -1.139e+01, -3.650e-02)
LowCD4 - AIDS      0.63465 (  9.115e-06,  4.419e+04)
LowCD4 - Death     0.01012 ( 8.238e-299, 1.242e+294)
AIDS - AIDS       -0.34837 ( -3.805e+40, -3.189e-42)
AIDS - Death       0.34837 (  3.189e-42,  3.805e+40)

-2 * log-likelihood:  1572.2

4. Construct the Rate Matrix

Rmsm <- msm::qmatrix.msm(hiv.msm, ci = "none")
round(Rmsm,3)

        Healthy LowCD4   AIDS Death
Healthy  -0.327  0.327  0.000 0.000
LowCD4    0.000 -0.645  0.635 0.010
AIDS      0.000  0.000 -0.348 0.348
Death     0.000  0.000  0.000 0.000

5. Embed the Transition Probability Matrix

Multistate-Model:

round(expm(Rmsm),3)

        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.525 0.388 0.088
AIDS      0.000  0.000 0.706 0.294
Death     0.000  0.000 0.000 1.000

Original Model Matrix:

round(mP,3)

        Healthy LowCD4  AIDS Death
Healthy   0.721  0.202 0.067 0.010
LowCD4    0.000  0.581 0.407 0.012
AIDS      0.000  0.000 0.750 0.250
Death     0.000  0.000 0.000 1.000

Markov Trace Comparison

# Original Model
c(1000, 0, 0, 0) %*% mP

     Healthy LowCD4 AIDS Death
[1,]     721    202   67    10

# MSM-Based Model
c(1000, 0, 0, 0) %*% pmatrix.msm(hiv.msm, t=1)

     Healthy LowCD4 AIDS Death
[1,]     721    202   67    10

Resulting expectation after 1 cycle is nearly identical.
The msm method has bits down at 4th decimal. Round off made them identical.
Only a tiny bit of error (721 vs. 721.0005) makes log numerically unstable converting from probability to rate.

Key Takeaways

By imposing some constraints on the underlying transitions, we were able to yield a generator matrix that makes sense!
Critically, our multistate model-based generator matrix closely approximates the original transition probability matrix, but does not imply negative rates that bring people back from death.

Key Takeaways

Multistate-Model:

round(logm(mP),3)

       [,1]   [,2]   [,3]   [,4]
[1,] -0.327  0.311  0.002  0.013
[2,]  0.000 -0.543  0.615 -0.072
[3,]  0.000  0.000 -0.288  0.288
[4,]  0.000  0.000  0.000  0.000

Original Model Matrix:

round(Rmsm,3)

        Healthy LowCD4   AIDS Death
Healthy  -0.327  0.327  0.000 0.000
LowCD4    0.000 -0.645  0.635 0.010
AIDS      0.000  0.000 -0.348 0.348
Death     0.000  0.000  0.000 0.000

Key Takeaways

With a reasonable generator matrix defined, we can now augment the model as we see fit:
- Add additional health states.
- Add new strategies with evidence on efficiacy drawn from the literature (e.g., hazard rates).
- Add accumulators and transition states.
- Change the cycle length.
Once done with the above, just embed the new matrix.

2. Solving for PSA Distributions

A Common Issue

We have defined our underlying model (either from bottom up or through backwards conversion) but need to define PSA distributions for model parameters.
Model draws on literature-based parameters that are reported as means, standard deviations, interquartile ranges, 95% confidence intervals, etc.

A Common Issue

Straightforward to obtain base case values from literature.
But how can we define PSA distributions based on limited information?

PSA Distributions

Parameter Type	Distribution
Probability	beta
Rate	gamma
Utility weight	beta
Right skew (e.g., cost)	gamma, lognormal
Relative risks or hazard ratios	lognormal
Odds Ratio	logistic

Example

Cost parameter reported in literature has interquartile range of $300-$750.
What are the parameters of a PSA distribution such that the 25th percentile is $300 and the 75th percentile is $750?

Some Options

Analytic formulas for some distributions (normal, gamma, beta).
Formulas take as their inputs the two values ($x_1$, $x_2$) and their associated quantiles ($p_1$,$p_2$).
Formulas return the PSA distribution parameters that match up to these values.

Some Options

ParameterSolver implemented in Windows, Link
See Cook (2010) and Vasco Grilo’s blog post for more.
The next few slides provide you with nearly all the tools you’ll need, however.

Example: Cost PSA Distribution

Suppose the interquartile range for a key cost variable is [300,750]
We may have obtained this from the literature, from tabulating cost data, or from expert opinions.

Example: Cost PSA Distribution

Suppose the interquartile range for a key cost variable is [300,750]
We may have obtained this from the literature, from tabulating cost data, or from expert opinions.

x1 = 300
p1 = 0.25

x2 = 750
p2 = 0.75

Analytic Solution: Uniform PSA

For a uniform distribution with minimum $a$ and maximum $b$

\[ a = \frac{p_2x_1 - p_1x_2}{p_2-p_1} \] \[ b = \frac{(1-p_1)x_2-(1-p_2)x_1}{p_2-p_1} \]

Analytic Solution: Uniform PSA

a = ((p2*x1) - (p1 * x2)) / (p2 - p1)
a

[1] 75

b = ((1 - p1)*x2 - (1 - p2)*x1) / (p2 - p1)
b

[1] 975

Analytic Solution: Uniform PSA

qunif(0.25, min = 75, max = 975)

[1] 300

qunif(0.75, min = 75, max = 975)

[1] 750

plot(density(runif(n = 1e5, min = a, max = b)),main = "")

Analytic Solution: Normal PSA

\[ \sigma = \frac{x_2 - x_1}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

\[ \mu = \frac{x_1\Phi^{-1}(p_2)-x_2\Phi^{-1}(p_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

Analytic Solution: Normal PSA

mu = (qnorm(p2)*x1 - qnorm(p1)*x2) / (qnorm(p2)-qnorm(p1))
mu

[1] 525

sigma = (x2 - x1) / (qnorm(p2) - qnorm(p1))
sigma

[1] 333.59

Analytic Solution: Normal PSA

qnorm(0.25, mean = 525, sd = 333.59)

[1] 300

qnorm(0.75, mean = 525, sd = 333.59)

[1] 750

plot(density(rnorm(n = 1e5,  mean = mu, sd = sigma)),main = "")

Analytic Solution: Lognormal PSA

Just take log of $x_1$ and $x_2$

\[ \sigma = \frac{\ln(x_2) - \ln(x_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

\[ \mu = \frac{\ln(x_1)\Phi^{-1}(p_2)-\ln(x_2)\Phi^{-1}(p_1)}{\Phi^{-1}(p_2)-\Phi^{-1}(p_1)} \]

Analytic Solution: Lognormal PSA

mu = (qnorm(p2)*log(x1) - qnorm(p1)*log(x2)) / (qnorm(p2)-qnorm(p1))
mu

[1] 6.1619

sigma = (log(x2) - log(x1)) / (qnorm(p2) - qnorm(p1))
sigma

[1] 0.67925

Analytic Solution: Lognormal PSA

qlnorm(0.25, mean = 6.1619, sd = 0.67925)

[1] 299.99

qlnorm(0.75, mean = 6.1619, sd = 0.67925)

[1] 749.98

plot(density(rlnorm(n = 1e5,  mean = mu, sd = sigma)),main = "")

Analytic Solution: Gamma PSA

A bit more involved as it involves finding the root of a function.

Analytic Solution: Gamma PSA

A bit more involved as it involves finding the root of a function.

gamma_fn <- function(alpha) {
    x1*qgamma(p2,shape = alpha, scale =1) - x2 * qgamma(p1, shape = alpha, scale = 1)
}
curve(gamma_fn, xlim = c(1,10), col = "blue", lwd = 1.5, lty=2)
abline(a=0,b=0)

Analytic Solution: Gamma PSA

Root (i.e., point where the function crosses the zero line) seems to be between 2 and 4.

gamma_fn <- function(alpha) {
    x1*qgamma(p2,shape = alpha, scale =1) - x2 * qgamma(p1, shape = alpha, scale = 1)
}
curve(gamma_fn, xlim = c(1,10), col = "blue", lwd = 1.5, lty=2)
abline(a=0,b=0)

Analytic Solution: Gamma PSA

We’ll search in this range for our $\alpha$ value.

gamma_fn <- function(alpha) {
    x1*qgamma(p2,shape = alpha, scale =1) - x2 * qgamma(p1, shape = alpha, scale = 1)
}
curve(gamma_fn, xlim = c(1,10), col = "blue", lwd = 1.5, lty=2)
abline(a=0,b=0)

Analytic Solution: Gamma PSA

alpha_ <- uniroot(gamma_fn,c(2,4))$root
alpha_

[1] 2.4558

calc_beta <- function(x1,p1,alpha) {
    x1 / qgamma(p1,alpha,1)
}

beta_ <- calc_beta(x1 = x1,  p1 = p1, alpha = alpha_)
beta_

[1] 230.16

Analytic Solution: Gamma PSA

qgamma(0.25,shape = 2.4558, scale = 230.16)

[1] 300

qgamma(0.75,shape = 2.4558, scale = 230.16)

[1] 750

plot(density(rgamma(n = 1e4, shape = alpha_, scale = beta_)),main = "")

Alternative: Optimization

General solution that can be used to solve for parameters of many common PSA distributions (e.g., beta, gamma, etc.).

\[\min_{\vec{\theta} \in R}\,\, (F(x_1|\vec{\theta})-p_1)^2+(F(x_2|\vec{\theta})-p_2)^2\]

Alternative: Optimization

Requires cumulative distribution function (CDF), $F$.
See Hans W Borchers’ great slides for tips on optimization in R.

Some Advice

We have found that while analytically correct, this optimization formula is not numerically stable.
Transfinite scaling of the probabilities from CDFs $F(.)$ stabilizes optimization.

Some Advice

We have found that while analytically correct, this optimization formula is not numerically stable.
Transfinite scaling of the probabilities from CDFs $F(.)$ stabilizes optimization.
$g(F(.))$ vs. $F(.)$

\[\min_{\vec{\theta} \in R}\,\, (g(F(x_1|\vec{\theta}))-p_1)^2+(g(F(x_2|\vec{\theta}))-p_2)^2\]

Stable Optimization

For probabilities, this is the $\text{logit}(p)=\log \frac{p}{1-p}=\log p - \log (1-p)$.
Tweaking parameters to converge can still happen.
Example of theory versus practice.

Stable Optimization

For probabilities, this is the $\text{logit}(p)=\log \frac{p}{1-p}=\log p - \log (1-p)$.
Tweaking parameters to converge can still happen.
Example of theory versus practice.

\[\min_{\vec{\theta} \in R} \sum_{i \in 1,2}\,\, (\text{logit} (F(x_i|\vec{\theta}))-\text{logit} (p_i))^2\]

Transfinite Example, Part 1

# Transfinite scaling
Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

# Function to minimize
norm <- function(x1, p1, x2, p2, shape,rate)
    (Tf(x1, shape, rate)-(log(p1)-log(1-p1)) )^2 +
    (Tf(x2, shape, rate)-(log(p2)-log(1-p2)) )^2

# Bundle it all together into a single function.
fn <- function(x) norm(x1, p1, x2, p2, x[1], x[2])

Transfinite Example, Part 2

# Run general-purpose optimization on the function.
gamma_optim <- 
  optim(c(0.5, 0.1), # initial parameter guesses
    fn,
    gr = function(x) pracma::grad(fn, x),
    lower=c(-Inf, 1e-5),
    method = "L-BFGS-B",
    control=list(factr=1e-10, maxit=100))$par

# Note: gamma_optim$par[2] is 1/beta
gamma_optim

[1] 2.4558556 0.0043449

Transfinite Example

Optimization returns $\alpha$ and $1/\beta$ for the gamma distribution example.

# Analytic formula
alpha_

[1] 2.4558

# Optimization
gamma_optim[[1]]

[1] 2.4559

# Analytic formula
beta_

[1] 230.16

# Optimization
1/gamma_optim[[2]]

[1] 230.16

Optimization for Other Distributions

Just swap in the distribution function for the distribution you want to match to (e.g., pbeta rather than pgamma).
Also make sure the supplied parameter names match the distribution you’re aiming for.

Optimization for Other Distributions

Just swap in the distribution function for the distribution you want to match to (e.g., pbeta rather than pgamma).
Also make sure the supplied parameter names match the distribution you’re aiming for.

Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

Optimization for Other Distributions

Just swap in the distribution function for the distribution you want to match to (e.g., pbeta rather than pgamma).
Also make sure the supplied parameter names match the distribution you’re aiming for.

Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

Summary on Distribution Fitting

Use analytical formulas if they exist.
Optimize on $\text{logit}$ scale.
Again, avoid $\log$ on probabilities, use log=TRUE.
Plot results for visual check.

Copula PSA Sampling

Sensitivity Analysis for Decision Models

It is standard in CEA analyses to explore sensitivity of model outputs and outcomes to variation in model inputs.
Many decision problems are complicated by dependencies between model inputs (e.g., costs, utilities, etc.).
These correlations may have meaningful impacts on model results and sensitivity.

HIV Model

In the HIV model discussed earlier, patients accrue two types of costs:

Direct medical costs associated with each health state (Healthy, LowCD4, AIDS)
Community medical costs associated with each health state.

HIV Model

These costs are likely correlated with each other in meaningful ways.
Similarly, other model parameters (e.g., state utility payoffs, hazard rates for treatment strategies, etc.) may also be correlated.

HIV Model

How can we reflect this correlation in n-way sensitivity analyses, such as PSAs?

What You Were Taught

Nearly all PSA exercises draw from independent distributions.
To fully explore the space of correlated parameters, you need a lot of PSA draws!

What You Were Taught

Independent PSA draws of two model parameters
- $c_{direct} \sim \text{gamma}(2.5, 1/230)$
- $c_{community} \sim \text{gamma}(1.5, 1/150)$

x1 <- rgamma(n = 1e3, shape = 2.5, rate  = 1/230.16)
x2 <- rgamma(n = 1e3, shape = 1.5, rate = 1/150)

tibble(x1 = x1, x2 = x2) %>% 
  ggplot(aes(x = x2, y = x1)) +
  geom_point(colour = "red") + 
  hrbrthemes::theme_ipsum_pub(base_family = "Arial") + 
  labs(x = "Direct Medical Costs", y = "Community Costs")

What We’ll Teach You

$\text{cor}(c_{direct},c_{community}) = 0.8$

# 1. Define the correlation matrix
  rho <- 0.8
  sigma <- matrix(c(1,rho,rho,1),nrow = 2, byrow=TRUE)

# 2. Perform a cholesky decomposition
  A = chol(sigma)

# 3. Generate iid standard normal pseudo random variables
  tildeY0 <- rnorm(n = 1e3, mean = 0, sd = 1)
  tildeY1 <- rnorm(n = 1e3, mean = 0, sd = 1)
  tildeY <- cbind(tildeY0,tildeY1)

  # Collect them 
  tildeY <- tildeY %*% A

# Use the standard normal disribution function to return the quantiles
  U <- pnorm(tildeY)

# Use the inverse distribution function to return the values
  Y <- U 
  cx1 <- qgamma(U[,1],shape = 2.5, rate  = 1/230.16)
  cx2 <- qgamma(U[,2],shape = 1.5, rate = 1/150)

tibble(x1 = cx1, x2 = cx2) %>% 
  ggplot(aes(x = x2, y = x1)) +
  geom_point() + 
  hrbrthemes::theme_ipsum_pub(base_family = "Arial") + 
  labs(x = "Direct Medical Costs", y = "Community Costs")

What We’ll Teach You

Red points: independent PSA draws
Black points: correlated PSA draws

tibble(x1 = x1, x2 = x2) %>% 
  ggplot(aes(x = x2, y = x1)) +
  geom_point(alpha = 0.5, colour = "red") + 
  geom_point(data = tibble(x1 = cx1, x2 = cx2) , aes(x = cx2, y = cx1),alpha = 0.25) +
  hrbrthemes::theme_ipsum_pub(base_family = "Arial") + 
  labs(x = "Direct Medical Costs", y = "Community Costs")  + 
  stat_ellipse(
    data = MASS::mvrnorm(n=1e3, mu = c(1500,750),Sigma = matrix(c(40000,0,0,20000),byrow=TRUE,nrow=2)) %>% 
      data.frame() %>% 
      set_names(c("x1","x2")) , 
    lty = 3, colour = "black"
  ) + 
  annotate("text", x = 750, y = 1500, label = "Traditional PSA draws (red points)\n very sparse in this region.",fontface =2, color = "red")

Copula-Based PSA Sampling

We can easily draw correlated PSA samples via copulas.
Copulas allow us to sample from the joint distribution of correlated parameters.

Copulas

Copula: a multivariate cumulative distribution function with uniform marginals.
Intuition: the binding “glue” between random variables.

Sklar (1959)

Sklar (1959) demonstrates that the joint distribution of two random variables can be represented as

Sklar (1959)

Sklar (1959) demonstrates that the joint distribution of two random variables can be represented as

\[ \mathrm{F}_{X_{1}, X_{0}}\left(x_1, x_0\right)=C_{X_{1}, X_{0}}\left(\mathrm{~F}_{X_{1} }\left(x_1\right), \mathrm{F}_{X_{0} }\left(x_0\right)\right) \]

Copula-Based PSA Draws

We’ll draw on copulas so we can sample our PSA as the joint distribution of uncertain parameters.
You can always just set the correlations to zero if you want independent draws!

\[ F(x_1,...,x_d) = C \big ( F_1(x_1), ..., F_d(x_d) \big ) \]

How Do We Sample?

Recall that the copula specifies multivariate CDF with uniform marginals.
What if our PSA distributions are not uniform (e.g., costs $\sim$ gamma, lognormal, etc.)?
Easy! Use the inverse transform method.
Idea: randomly sample from quantiles in the PSA distribution.

How Do We Draw PSA Values

Requires drawing a random number between 0 and 1 (the quantile)
Requires an inverse cumulative distribution function $F^{-1}(.)$
You can do this easily in R or Excel!

Normal Distribution: PDF to CDF

library(tidyverse)
library(data.table)
library(gganimate)
set.seed(12)
df_norm <- tibble(x=rnorm(10000,mean = 10, sd=3))
df_lnorm <- tibble(x = rlnorm(10000,mean=log(3),sd=0.01)); plot(density(df_lnorm$x))
df_gamma <- tibble(x = rgamma(10000,shape = 44.4, scale=22.5)); plot(density(df_gamma$x))
df_unif <- tibble(x = runif(1e6,min = 0.3,max=0.4)); plot(density(df_unif$x))

animate_pdf_to_cdf <- function(df, scene_to_show = NULL, q_ = c(0.01,  0.25, 0.5, 0.75,  0.99), breaks_to_show = NULL ) {
  pdf <-
    {df %>% ggplot(aes(x = x)) + geom_density()} %>%
    ggplot_build() %>% pluck("data") %>% pluck(1) %>%
    as_tibble() %>%
    select(x,y) %>%
    data.table(key = "x")
  cdf <-
    {df %>% ggplot(aes(x = x)) + stat_ecdf()} %>%
    ggplot_build() %>% pluck("data") %>% pluck(1) %>%
    as_tibble() %>%
    select(x,y) %>%
    data.table(key = "y")

  
  q <- q_ %>%
    map_dbl(~{
      quantile(df$x,.x)
    }) %>%
    set_names(q_) %>%
    data.frame() %>%
    set_names("x_") %>%
    rownames_to_column(var = "quantile") %>%
    as.data.table(key = "x_")

  scene2 <-
    pdf[q,roll="nearest"] %>%
    as_tibble() %>%
    select(y,quantile) %>%
    right_join(pdf,c("y")) %>%
    arrange(x,y) %>%
    mutate(scene = 2) %>%
    mutate(quantile = glue::glue("{as.numeric(quantile)}"))

  q2 <-
    q %>%
    as_tibble() %>%
    mutate(quantile = as.numeric(quantile)) %>%
    data.table(key="quantile")

  scene3 <-
    cdf[q2,roll="nearest"] %>%
    as_tibble() %>%
    set_names(c("x","quantile","x_")) %>%
    select(-x_) %>%
    right_join(cdf,"x") %>%
    arrange(x,y)  %>%
    mutate(scene = 3) %>%
    mutate(quantile = glue::glue("{as.numeric(quantile)}")) %>%
    ungroup() %>%
    filter(!is.na(as.numeric(quantile)) | row_number() %in% sample(1:nrow(.),1000)); dim(scene2)

  p_df <-
    scene2 %>%
    mutate(scene = 1) %>%
    bind_rows(scene2) %>%
    bind_rows(scene3) %>%
    bind_rows(scene3 %>% mutate(scene=4)) %>%
    #bind_rows(scene3 %>% mutate(scene=5)) %>%
    #bind_rows(scene3 %>% mutate(scene=6)) %>%
    mutate(quantile = ifelse(quantile=="NA",NA, quantile))

  p_df_label =
    p_df %>% filter(!is.na(quantile)) %>%
    mutate(y = ifelse(scene==1,0,y)) %>%
    mutate(x = ifelse(scene==4,min(scene3$x),x)) %>%
    mutate(y = ifelse(scene==6,0,y)) %>%
    mutate(quantile = glue::glue(" q{quantile}"))

  if (!is.null(scene_to_show)) {
    p_df <- 
      p_df %>% filter(scene==scene_to_show)
    p_df_label <- 
      p_df_label %>% filter(scene==scene_to_show)
    p <-
      p_df %>% ggplot(aes(x = x, y = y))+ geom_line() + geom_text(data = p_df_label , aes(label = quantile, group= quantile),size=3)+
      scale_y_continuous(breaks = breaks_to_show) +
      scale_x_continuous(breaks = q$x_, labels = round(q$x_,2),guide = guide_axis(n.dodge=3)) +
      ggthemes::theme_clean(base_size=15) +
      labs(x = "Parameter Value", y = "") 
  } else {
      p <-
        p_df %>% ggplot(aes(x = x, y = y))+ geom_line() + geom_text(data = p_df_label , aes(label = quantile, group= quantile))+
        scale_y_continuous(breaks = seq(0,1,0.1)) +
        scale_x_continuous(breaks = q$x_, labels = round(q$x_,2)) +
        ggthemes::theme_clean() +
        labs(x = "Parameter Value", y = "") + #facet_wrap(~scene, scales="free"); p
        gganimate::transition_states(scene) + gganimate::view_follow()
  }

  return(p)

}

df_norm %>% animate_pdf_to_cdf(scene = 2)
ggsave(here::here("slides/images/pdf-to-cdf-norm.png"),height=4, width=5)
df_norm %>% animate_pdf_to_cdf(scene = 3,seq(0,1,0.1)) 
ggsave(here::here("slides/images/pdf-to-cdf-norm4.png"),height=4, width=5)
set.seed(23)
df_norm %>% animate_pdf_to_cdf(scene = 3, q = round(runif(5, min = 0, max = 1),3))
ggsave(here::here("slides/images/pdf-to-cdf-norm_random1.png"),height=4, width=5)
anim_norm <- df_norm %>% animate_pdf_to_cdf(breaks_to_show = seq(0,1,0.1))
anim_save(here::here("slides/images/pdf-to-cdf-norm.gif"), anim_norm)

anim_gamma1 <- df_gamma %>% animate_pdf_to_cdf(scene = 2)
ggsave(here::here("slides/images/pdf-to-cdf-gamma.png"),height=5, width=5)
anim_gamma3 <- df_gamma %>% animate_pdf_to_cdf(scene = 4)
ggsave(here::here("slides/images/pdf-to-cdf-gamma4.png"),height=5, width=5)
anim_gamma <- df_gamma %>% animate_pdf_to_cdf()
anim_save(here::here("slides/images/pdf-to-cdf-gamma.gif"), anim_gamma)

anim_lnorm1 <- df_lnorm %>% animate_pdf_to_cdf(scene = 2)
ggsave(here::here("slides/images/pdf-to-cdf-lnorm.png"),height=5, width=5)
anim_lnorm3 <- df_lnorm %>% animate_pdf_to_cdf(scene = 4)
ggsave(here::here("slides/images/pdf-to-cdf-lnorm4.png"),height=5, width=5)
anim_lnorm <- df_lnorm %>% animate_pdf_to_cdf()
anim_save(here::here("slides/images/pdf-to-cdf-lnorm.gif"), anim_lnorm)

set.seed(123)
r_HS =  tibble(x = rgamma(1e4,shape = 30,rate = 200),param="r_HS")
r_HD    = tibble(x = rgamma(1e4,shape = 60, rate = 10000),param = "r_HD")
hr_S =  tibble(x = rlnorm(1e4,log(3), sdlog = 0.01),param = "hr_S")
hr_HS_trtB  = tibble(x=rlnorm(1e4,log(0.96), sdlog = 0.02),param = "hr_HS_trtB")
hr_HS_trtC  = tibble(x=rlnorm(1e4,log(0.9), sdlog = 0.02),param = "hr_HS_trtC")
hr_HS_trtD =    tibble(x=rlnorm(1e4,log(0.92), sdlog = 0.02),param = "hr_HS_trtD")
hr_HS_trtE =    tibble(x=rlnorm(1e4,log(0.92), sdlog = 0.02), param = "hr_HS_trtE")
u_H = tibble(x=rbeta(1e4,shape1 = 200, shape2 = 3), param = "u_H")
u_S = tibble(x=rbeta(1e4,shape1 = 130, shape2 = 45), param = "u_S")
c_S  = tibble(x=rgamma(1e4,shape = 44.4, scale = 22.5), param = "c_S")
c_trtA  =tibble(x= rgamma(1e4,shape = 12.5, scale = 2), param = "c_trtA")
c_trtB  =tibble(x= rgamma(1e4,shape = 12, scale = 83.3), param = "c_trtB")
c_trtC  =tibble(x= rgamma(1e4,shape = 36.144, scale = 83), param = "c_trtC")
c_trtD  =tibble(x= rgamma(1e4,shape = 14.458, scale = 83), param = "c_trtD")
c_trtE  =tibble(x= rgamma(1e4,shape = 60.24, scale = 83), param = "c_trtE")


params_lut <- c("r_HS = gamma(shape = 30,rate = 200)", "r_HD = gamma(shape = 60, rate = 10000)", "hr_S = lognormal(meanlog = log(3), sdlog = 0.01)", "hr_HS_trtB = lognormal(meanlog = log(0.96), sdlog = 0.02)", "hr_HS_trtC = lognormal(meanlog = log(0.88), sdlog = 0.02)", "hr_HS_trtD = lognormal(meanlog = log(0.92), sdlog = 0.02)", "hr_HS_trtE = lognormal(meanlog = log(0.92), sdlog = 0.02)",
               "u_H = beta(shape1 = 200, shape2 = 3)","u_S = beta(shape1 = 130, shape2 = 45)","c_S = gamma(shape = 44.4, scale = 22.5)", "c_trtA = gamma(shape = 12.5, scale = 2)","c_trtB = gamma(shape = 12, scale = 83.3)","c_trtC = gamma(shape = 36.144, scale = 83)","c_trtD = gamma(shape = 18.67, scale = 83)","c_trtE = gamma(shape = 60.24, scale = 83)")


params <- list(r_HS, r_HD, hr_S, hr_HS_trtB, hr_HS_trtC, hr_HS_trtD, hr_HS_trtE,
               u_H,u_S,c_S, c_trtA,c_trtB,c_trtC,c_trtD,c_trtE)

pdf_params <- 
  params %>% 
  map2(.,params_lut,~({
    .x %>% animate_pdf_to_cdf(scene = 1) + ggtitle(.y)
    ggsave(here::here(glue::glue("slides/images/{.x$param[1]}.png")),height=4, width=5)
  })) %>% 
  set_names(c("r_HS:gamma(shape = 30,rate = 200)", "r_HD=gamma(shape = 60, rate = 10000)", "hr_S=lognormal(meanlog = log(3), sdlog = 0.01)", "hr_HS_trtB=lognormal(meanlog = log(0.96), sdlog = 0.02)", "hr_HS_trtC=lognormal(meanlog = log(0.88), sdlog = 0.02)", "hr_HS_trtD=lognormal(meanlog = log(0.92), sdlog = 0.02)", "hr_HS_trtE=lognormal(meanlog = log(0.92), sdlog = 0.02)",
               "u_H=beta(shape1 = 200, shape2 = 3)","u_S=beta(shape1 = 130, shape2 = 45)","c_S=gamma(shape = 44.4, scale = 22.5)", "c_trtA=gamma(shape = 12.5, scale = 2)","c_trtB=gamma(shape = 12, scale = 83.3)","c_trtC=gamma(shape = 36.144, scale = 83)","c_trtD=gamma(shape = 18.67, scale = 83)","c_trtE=gamma(shape = 60.24, scale = 83)"))

cdf_params <- 
  params %>% 
  map2(.,params_lut,~({
    set.seed(23)
    .x %>% animate_pdf_to_cdf(scene = 3,q_ = round(runif(5, min = 0, max = 1),3),breaks_to_show = seq(0,1,0.1)) + ggtitle(.y)
    ggsave(here::here(glue::glue("slides/images/cdf{.x$param[1]}.png")),height=4, width=5)
  })) %>% 
  set_names(c("r_HS", "r_HD", "hr_S", "hr_HS_trtB", "hr_HS_trtC", "hr_HS_trtD", "hr_HS_trtE",
               "u_H","u_S","c_S", "c_trtA","c_trtB","c_trtC","c_trtD","c_trtE"))

anim_norm

Normal Distribution: CDF

Drawing from an Arbitrary PSA Distribution

Define PSA distribution parameters (e.g., $\mu,\sigma$ if normal, etc.)
Draw a uniform random number between 0 and 1.
Feed this number through the quantile function for the specified distribution (e.g., qnorm(), qgamma(), etc.). This returns the value that maps to the randomly drawn quantile.
Do this for as many PSA samples as you need!

Copula-Based Sampling

Sampling from a copula adds one more step to this process.
We must also define a correlation matrix Sigma ($\Sigma$).
We then draw uniform random numbers from a multivariate uniform distribution.
Once drawn, we can feed each PSA variable back through its repective inverse CDF to obtain the final PSA sample.

Copula-Based Sampling

Define PSA distributions for each of $n$ uncertain parameters.
Define an $n \times n$ correlation matrix $\Sigma$ for the joint distribution of the PSA parameters.
Draw an $n$-column multivariate uniform sample.
Apply the inverse transform method (as necessary) to each column.
Do this for as many PSA samples as you need!

Putting it All Together

IQR Community Care Costs: 750 to 1250
- cost_cc $\sim$ lognormal(mu_cc, sigma_cc).
IQR Direct Medical Care Costs: 1,000 - 3,000
- cost_dc $\sim$ gamma(shape_dc, rate_dc)
COV(Community cost, Direct cost) = 0.8

1. Define PSA distributions

IQR Community Care Costs: 750 to 1250
- cost_cc $\sim$ lognormal(mu_cc, sigma_cc).

mu_cc = (qnorm(0.75)*log(750) - qnorm(0.25)*log(1250)) / (qnorm(0.75)-qnorm(.25)); mu_cc

[1] 6.8755

sigma_cc = (log(1250) - log(750)) / (qnorm(.75) - qnorm(.25)); sigma_cc

[1] 0.37868

1. Define PSA distributions

IQR Direct Medical Care Costs: 1,000 - 3,000
- cost_dc $\sim$ gamma(shape_dc, rate_dc)

x1 = 1000
x2 = 3000
p1 = 0.25
p2 = 0.75

# Transfinite scaling
Tf <- function(x, shape, rate)
    pgamma(x,shape,rate,log=TRUE) - 
    pgamma(x,shape,rate,log=TRUE, lower.tail = FALSE)

# Function to minimize
norm <- function(x1, p1, x2, p2, shape,rate)
    (Tf(x1, shape, rate)-(log(p1)-log(1-p1)) )^2 +
    (Tf(x2, shape, rate)-(log(p2)-log(1-p2)) )^2

# Bundle it all together into a single function.
fn <- function(x) norm(x1, p1, x2, p2, x[1], x[2])

# Run general-purpose optimization on the function.
gamma_optim <- 
  optim(c(0.5, 0.1), # initial parameter guesses
    fn,
    gr = function(x) pracma::grad(fn, x),
    lower=c(-Inf, 1e-5),
    method = "L-BFGS-B",
    control=list(factr=1e-10, maxit=100))$par

shape_dc = gamma_optim[[1]]; shape_dc

[1] 1.7904

rate_dc = gamma_optim[[2]]; rate_dc

[1] 0.00080927

2. Define correlation matrix $\Sigma$

rho_cc_dc = 0.8
sigma <- matrix(c(1,rho_cc_dc,rho_cc_dc,1),nrow = 2, byrow=TRUE,
                dimnames=list(c("cost_cc","cost_dc"),c("cost_cc","cost_dc")))
sigma

        cost_cc cost_dc
cost_cc     1.0     0.8
cost_dc     0.8     1.0

3. Draw multivariate uniform sample.

To get multivariate uniform, first draw from multivariate normal (MASS::mvrnorm()).
Then, apply the values to the normal CDF (pnorm()) to get the quantiles.

3. Draw multivariate uniform sample.

To get multivariate uniform, first draw from multivariate normal (MASS::mvrnorm()).
Then, apply the values to the normal CDF (pnorm()) to get the quantiles.

# Define PSA Sample Size
PSA_sample_size = 1e6
U <- pnorm(MASS::mvrnorm(n = PSA_sample_size, Sigma = sigma, mu = c(0,0)))

# Alternative Version: Cholesky decomposition

#  Perform a cholesky decomposition
  A = chol(sigma)

#  Generate iid standard normal pseudo random variables
  tildeY0 <- rnorm(n = PSA_sample_size, mean = 0, sd = 1)
  tildeY1 <- rnorm(n = PSA_sample_size, mean = 0, sd = 1)
  tildeY <- cbind(tildeY0,tildeY1)

  # Collect them 
  tildeY <- tildeY %*% A

# Use the standard normal disribution function to return the quantiles
  U <- pnorm(tildeY)

4. Inverse Transform

# Inverse transform using parameters solved for earlier
cost_cc <- qlnorm(U[,1], mu_cc, sigma_cc)
cost_dc <- qgamma(U[,2], shape_dc, rate_dc)

# Make sure the marginal distributions check out
quantile(cost_cc,c(0.25,0.75))

    25%     75% 
 750.27 1249.99

quantile(cost_dc,c(0.25,0.75))

    25%     75% 
 999.21 2997.30

# Check that correlation is there
cor(cost_cc,cost_dc)

[1] 0.78332

Your Turn!

References

Chhatwal, Jagpreet, Suren Jayasuriya, and Elamin H. Elbasha. 2016. “Changing Cycle Lengths in State-Transition Models: Challenges and Solutions.” Medical Decision Making 36 (8): 952–64. https://doi.org/10.1177/0272989X16656165.

Cook, John D. 2010. “Determining Distribution Parameters from Quantiles.” Working paper.

Higham, Nicholas J. 2008. Functions of Matrices: Theory and Computation. SIAM.